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We study the effects of quenched disorder on the first-order phase transition in the two-dimensional 
three-color Ashkin-Teller model by means of large-scale Monte Carlo simulations. We demonstrate 
that the first-order phase transition is rounded by the disorder and turns into a continuous one. Us¬ 
ing a careful finite-size-scaling analysis, we provide strong evidence for the emerging critical behavior 
of the disordered Ashkin-Teller model to be in the clean two-dimensional Ising universality class, 
accompanied by universal logarithmic corrections. This agrees with perturbative renormalization- 
group predictions by Cardy. As a byproduct, we also provide support for the strong-universality 
scenario for the critical behavior of the two-dimensional disordered Ising model. We discuss conse¬ 
quences of our results for the classification of disordered phase transitions as well as generalizations 
to other systems. 
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I. INTRODUCTION 

The Imry-Ma criterion^ is one of the key results on 
phase transitions in disordered systems. It governs 
the stability of macroscopic phase coexistence against 
quenched random disorder that locally favors one phase 
over the other. By comparing the energy gain due to 
the disorder with the energy cost of a domain wall, Imry 
and Ma showed that disorder destroys phase coexistence 
by domain formation in dimensions d < As a con¬ 
sequence, infinitesimal disorder rounds first-order phase 
transitions in d < 2 as Aizenman and Wehr— later proved 
rigorously as a theorem.— 

These results raise the important question of what is 
the fate of a first-order transition that is destroyed by 
disorder? Is it a continuous transition, does an interme¬ 
diate phase appear, or is the sharp transition, perhaps, 
completely destroyed via smearing? If the transition be¬ 
comes continuous, what is the critical behavior? Is it 
accompanied by pretransitional singularities due to rare 
regions, as is the case at “generic” critical points in disor¬ 
dered systems (see, e.g., Refs. H andH)? These questions 
have recently reattracted considerable attention, in par¬ 
ticular in the context of zero-temperature quantum phase 

transitions^Ai 

It turns out, however, that these questions remain un¬ 
resolved even for a simple prototypical classical phase 
transition, viz., the transition in the two-dimensional fer¬ 
romagnetic Ashkin-Teller model. The TV-color Ashkin- 
Teller modeli^rJ^ consists of N Ising models, coupled via 
their energy densities. In the absence of disorder and 
for N > 2, this system features a first-order phase tran¬ 
sition between a paramagnetic high-temperature phase 
and a ferromagnetic (Baxter) phase a low temperatures. 
According to the Imry-Ma criterion, or equivalently the 


Aizenman-Wehr theorem, this first-order transition can¬ 
not survive the introduction of weak disorder in the form 
of random bonds or bond or site dilution. Murthyi^ and 
Cardyii analyzed this problem by means of perturba¬ 
tive renormalization group calculations which predicted 
that the first-order transition is rounded to a continuous 
transition in the universality class of the two-dimensional 
clean Ising model, apart from logarithmic corrections. 
However, recent numerical simulations of a random-bond 
three-color Ashkin-Teller mode l^^d^ disagreed with these 
predictions. They found nonuniversal critical exponents 
that vary with disorder strength and differ from the clean 
Ising exponents. Moreover, the reported value of the cor¬ 
relation length exponent v violates the inequality dv > 2 
due to Chayes et al<^ 

To resolve these contradicting results, we perform 
large-scale high-accuracy Monte Carlo simulations of the 
two-dimensional three-color Ashkin-Teller model. We 
consider two types of quenched disorder, random bonds 
as well as site dilution. Our data provide strong evidence 
that the emerging critical behavior is universal and in the 
clean Ising universality class, as predicted by the renor¬ 
malization group calculations]It is also accompanied 
by logarithmic corrections analogous to those found in 
the disordered two-dimensional Ising model. 

The rest of our paper is organized as follows: In Sec. 
El we introduce the A^-color Ashkin-Teller model and 
discuss its properties in the absence of disorder. We 
then briefly summarize the results of Cardy’s renormal¬ 
ization group theory. In Sec. uni we explain our Monte 
Carlo method, and we give an overview over the simu¬ 
lation parameters. Section m is devoted to the numer¬ 
ical results for the clean, site-diluted, and random-bond 
Ashkin-Teller models. As a byproduct, our data provide 
additional support for the strong-universality scenario for 
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the two-dimensional disordered Ising model. We con¬ 
clude in Sec.|V]by discussing consequences of our results 
for the classification of disordered phase transitions as 
well as generalizations to other systems. 

II. MODEL AND THEORY 

The two-dimensional A^-color Ashkin-Teller modeU^^— 
is a generalization of the original model proposed by 
Ashkin and Teller— (which corresponds to the N = 2 
case). It consists of N identical Ising models, coupled 
via their energy densities. The Hamiltonian of the clean 
model reads 

N 

a=l (ij) a<fi (ij) 

Here, i and j denote the sites of a regular square lattice of 
sites, and the corresponding sum is over pairs of near¬ 
est neighbors, a is the “color” index that distinguishes 
the N Ising models, and Sf = ±1 are the usual classical 
Ising variables. We are interested in the regime in which 
both the Ising interaction J and the four-spin interaction 
K are positive. The strength of the coupling between the 
Ising models can be parameterized by the dimensionless 
ratio e = K/J. Note that the Hamiltonian ([1]) is self 
dual for the case of iV = 2 colors; and this property has 
been used to find the exact location of the phase transi¬ 
tion in the clean and disordered models.— For N > 2, 
the Hamiltonian ([T]) is not self-dual. Self-duality can be 
restored, however, by including higher-order terms with 
up to 2N spins— We have not done this in our work, 
mainly to keep our results quantitatively comparable to 
other simulation o^^dSiiQ the literature. 

The properties of the clean Ashkin-Teller model have 
been studied in great detail. The two-color model {N = 
2) features a continuous transition with nonuniversal, 
continuously varying exponents between a paramagnetic 
high-temperature phase and an ordered phase at low tem¬ 
peratures (see, e.g.. Ref. [13 and references therein). In 
contrast, this transition is of first order for TV >2, which 
is the case we are interested in.— 

Quenched disorder can be introduced into the Hamil¬ 
tonian o in several ways. We consider both site dilution 
and bond randomness. In the former case, a fraction p 
of the lattice sites is removed at random (the 5”“ for all 
colors a are removed at such vacancy sites). The inter¬ 
actions between the remaining sites retain their uniform 
values J and K. In the case of bond randomness, the 
Ising couplings Jij between neighboring sites i and j be¬ 
come independent random variables drawn from some 
probability distribution P{J) which we take to be a bi¬ 
nary distribution 

PiJ) = cSiJ - A) + {1 - c)6iJ - Ji) (2) 

with Jfi > Ji > 0. Here, c is the concentration of the 
stronger bonds. The four-spin couplings Kij are either 



coupling strength e 


FIG. 1. (Color online) Cardy’s renormalization group trajec¬ 
tories on the critical surface (coupling strength e vs. disor¬ 
der strength A) for A = 3. The trajectories initially flow 
towards the first-order region at strong coupling. However, 
they eventually curl back towards the clean Ising fixed point 
at e = A = 0. 

taken to be uniform or they are slaved to the Ising inter¬ 
actions on the same bond via Kij = eJij with constant 
e. Both site dilution and random bonds are realizations 
of random-Tc disorder, i.e., disorder that does not break 
any of the spin symmetries but changes the local ten¬ 
dency towards the high-temperature or low-temperature 
phases. Thus, if the system undergoes a continuous phase 
transition, both types of disorder should lead to the same 
universality class. 

Murthyi^ and Cardyii applied a perturbative renor¬ 
malization group to a continuum version of the two- 
dimensional A-color Ashkin-Teller model. This analysis 
benefits from the fact that the first-order phase transition 
in the clean model is fluctuation-driven. As a result, the 
renormalization group is controlled in the limit of small 
inter-color coupling and weak disorder. Cardy found the 
renormalization group trajectories on the critical surface 
in closed form. In terms of the coupling strength e and 
the dimensionless disorder strength A, they read 

e = const X (A/e)^^-^)/^ exp(-2A/TVe) (3) 

A few characteristic trajectories are shown in Fig. [TJ For 
weak bare (initial) disorder, the trajectories first run to¬ 
wards the strong-coupling region e 1 where the tran¬ 
sition would turn first order. However, they eventually 
turn around and curl back towards the clean Ising fixed 
point at e = A = 0. This not only implies that the tran¬ 
sitions has become continuous, in agreement with the 
Imry-Ma criterion, it also means that the critical behav¬ 
ior is in the clean Ising universality class. A more detailed 
analysis of the renormalization group equations produces 
additional logarithmic corrections to the leading Ising 
power laws, similar to those found in renormalization 
group approaches to the disordered Ising model— 
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Furthermore, the large excursions of the renormaliza¬ 
tion group trajectories for small bare disorder strength 
imply a very slow crossover from the first-order tran¬ 
sition of the clean Ashkin-Teller model to the critical 
point of the disordered system. This crossover is espe¬ 
cially interesting because d = 2 is the marginal dimen¬ 
sionality for the Aizenman-Wehr theorem. (First-order 
transitions are destroyed by randomness for d < 2 while 
they can survive for d > 2.) If the clean system has 
a strong hrst-order transition, the breakup length Li,, 
beyond which randomness becomes important increases 
very rapidly with decreasing disorder strength A. In fact, 
for weak disorder, it is expecte d^^’^° to follow the expo¬ 
nential Lb ~ exp(const/A^). This implies that enormous 
system sizes are necessary to reach the asymptotic regime 
if the clean first-order transition is strong and the disor¬ 
der is weak. 


III. MONTE CARLO SIMULATIONS 

A. Overview 

To resolve the discrepancy between the renormaliza¬ 
tion group predictions outlined above and the recent nu¬ 
merical results of Refs. [H and [l^, we perform large- 
scale high-accuracy Monte Carlo simulations of two- 
dimensional three-color Ashkin-Teller models with site 
dilution and/or bond randomness. 

As we are interested in the critical behavior, a cluster 
algorithm is required to reduce the critical slowing down 
close to the phase transition. We employ a Wolff em¬ 
bedding algorithm similar to that used by Wiseman and 
DomanyS^ for the two-color Ashkin-Teller model. Its ba¬ 
sic idea is simple. Imagine fixing all and spins. 
Then, the Hamiltonian © is equivalent to an (embed¬ 
ded) Ising model for the spins with effective interac¬ 
tions 

Jff = J + eJ Sf + . (4) 

Simulating this Ising model using any valid Monte Carlo 
algorithm establishes detailed balance between all states 

(2) (S’) 

with the same fixed and S') . We can construct 

and simulate analogous embedded Ising models to up- 
(2) f3l 

date S') ' and S) '. By combing Monte Carlo updates for 
all three embedded Ising models we arrive at a valid algo¬ 
rithm (fulfilling ergodicity and detailed balance between 
all states) for the entire Hamiltonian ([1]). 

To simulate the embedded Ising models, we use the ef¬ 
ficient Wolff and Swendsen-Wang cluster algorithms 
They are only valid if all interactions are ferromagnetic, 
i.e., if all Jff > 0. This is fulfilled as long as the coupling 
strength |ef < 1/{N — 1). In our case of three colors, e 
therefore must not exceed 1/2. 

Finding the averages, variances, and distributions of 
observables in disordered systems requires the simulation 


of many samples with different disorder configurations. 
For optimal performance, one must therefore carefully 
choose the number Ug of samples (i.e., disorder configu¬ 
rations) and the number Um of measurements during the 
simulation of each samplei^^— Assuming statistical in¬ 
dependence between measurements (quite possible with 
a cluster algorithm), the total variance cr/ of a particular 
observable (thermodynamically and disorder averaged) 
can be estimated as 

O'* = + <^'Lln^)lns (5) 

where cr^ is the disorder-induced variance between sam¬ 
ples and is the variance of measurements within each 
sample. As the numerical effort is roughly proportional 
to rimns (neglecting equilibration for the moment), it is 
clear that the best value of rim is quite small. One might 
even be tempted to measure only once per sample. How¬ 
ever, with too few measurements, the majority of the 
computer time would be spent on equilibration. These 
requirements can be balanced by using large numbers ris 
of disorder configurations (ranging from several 10000 to 
several million in our case) and rather short runs with 
a few hundred Monte Carlo measurements per sample. 
Note that such short runs lead to biases in several ob¬ 
servables, at least if the usual estimators are employed. 
These biases can be corrected by improved estimators as 
is discussed in Appendix 0 

Based on these ideas, we develop two independent 
Monte Carlo codes, one (referred to as code A) mainly 
employed for simulating the site-diluted Ashkin-Teller 
model, and the other one (code B) used for the random- 
bond case. 


B. Site-diluted simulations 

All site-diluted simulations use code A. We study im¬ 
purity concentrations p = 0 (the clean case), 0.05, 0.1, 
0.2, and 0.3. For comparison, the lattice percolation 
threshold is at Pc = 0.407253. The Ising interaction J is 
fixed at unity while the coupling strength e = K/ J takes 
values 0 (the Ising limit), 0.05, 0.1, 0.2, 0.3, and 0.5. 
The system is tuned through the transition by changing 
the temperature T. Lattice sizes range from 25^ sites to 
1600^ sites (2240^ sites for the Ising case, e = 0) with pe¬ 
riodic boundary conditions. Data are averaged over up 
to 4 million disorder configurations for the smaller sys¬ 
tems and over up to 500,000 configurations for the largest 
ones. This leads to small statistical errors of the data. 

For site-diluted systems, we combine the Wolff single¬ 
cluster updates with Swendsen-Wang multi-cluster up¬ 
dates to equilibrate small isolated clusters of lattice sites 
that can occur for larger dilutions. Specifically, a full 
Monte Carlo sweep consists of a Swendsen-Wang sweep 
(for each color) followed by a Wolff sweep. (A Wolff sweep 
is defined as a number of cluster flips such that the total 
number of flipped spins per color is equal to the num¬ 
ber of sites.) To verify our codes we have also compared 
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the results to those of conventional Metropolis single-spin 
updates 1 ^ 

To estimate the equilibration times, we compare runs 
with “hot start” (initial spin values are completely ran¬ 
dom) and “cold start” (initially, all spins S'" = 1). Char¬ 
acteristic equilibration times range from less than 10 
sweeps for linear system size L = 50 to about 40 sweeps 
for system size L = 1600. In our production runs, we 
therefore employ equilibration periods of 60 to 100 sweeps 
and measurement periods of another 100 to 200 sweeps, 
with measurements taken after every sweep. Using these 
parameters, the results of runs with hot and cold starts 
agree within our small statistical errors. 

Note that simulations of the clean Ashkin-Teller model 
close to its strong first-order phase transition require 
longer equilibration times to overcome the supercritical 
slowing down associated with first-order transitions. De¬ 
tails will be given in Sec. IIV Al 


C. Random-bond simulations 


Using code A, we study the random-bond Ashkin- 
Teller model with the binary bond distribution ([2]). The 
Ising interactions take the values A = 2 or J; = 0.5, 
each with a probability of 0.5. The four-spin interactions 
are slaved to the Ising interactions via Kij = eJij with 
constant coupling strength e. We explore the cases e = 0 
(random-bond Ising model), 0.1, 0.2, and 0.5. Lattice 
sizes range from 35^ to 1120^ sites with periodic bound¬ 
ary conditions. The numbers of disorder configurations 
for each parameter set range from 10^ for the largest sys¬ 
tems to 10® for the smallest ones. 

Otherwise, the random-bond simulations are analo¬ 
gous to the site-diluted ones: Each full Monte Carlo 
sweep is a combination of a Wolff sweep and a Swendsen- 
Wang sweep. Each system is equilibrated from a hot start 
using 100 full sweeps, the measurement period is another 
100 sweeps. 

In addition, we use code B to perform simulations of 
a random-bond Ashkin-Teller model with uniform, non- 
random four-spin interaction K. Specifically, the Ising 
interactions take the values Jh = 6/5 and J/ = 4/5, each 
chosen with a probability of 0.5. This implies that the av¬ 
erage Ising interaction is J = {Jh + Ji)/2 = 1. The four- 
spin interactions are non-random and given hy K = eJ. 
Because the effective interactions O appearing in the 
embedded Wolff algorithm must be positive for both val¬ 
ues of the Ising interaction, the coupling strength is re¬ 
stricted to e < 2/5. In our simulations, e will be fixed at 
0.1. We simulate lattice sizes ranging from 24^ to 1600^ 
sites with periodic boundary conditions. The numbers 
of disorder configurations range from 10'^ for the largest 
system to 10® for the smallest one. Each system is equili¬ 
brated from a cold start using 200 full Wolff sweeps, the 
measurement period is also 200 Wolff sweeps per temper¬ 
ature; and we take measurements after every four sweeps. 


D. Observables 


During the simulations, we calculate various thermo¬ 
dynamic quantities such as the energy E = [(e)]dis and 
the magnetization M = [(TO)]dis- Here e and m stand 
for individual energy and magnetization measurements, 
and (...) is the canonical thermodynamic average (which 
is approximated by the Monte Carlo average over 
measurements). The average [. ..]dis over the disor¬ 
der distribution is approximated by the average over Ug 
samples. Specific heat and magnetic susceptibility are 
calculated from the fluctuations of e and m as C = 
(LVT2)[(e2) - (e)2]di. and y = {LyT)[{m^) - {m)^U. 
We also measure the product order parameter (or “po¬ 
larization”) Mp = [(mp)]dis with rup = (1/L^) Si •S'fS'f 
for two different colors a and /3. The corresponding sus¬ 
ceptibility reads Xp = {L'^/T)[{ml) - (mp)2]dis- 

Magnetization and susceptibility are averaged over the 
three colors for increased accuracy, and all quantities are 
normalized “per spin”. Analogously, the product order 
parameter and its susceptibility are averaged over the 
three possible pairs of colors. The statistical errors of 
all thermodynamic quantities are estimated from their 
fluctuations between disorder configurations. 

In addition, we calculate several quantities whose scale 
dimension is zero which makes them particularly suitable 
for a finite-size scaling analysis. The first such quantity 
is the Binder cumulant of the magnetization. In a disor¬ 
dered system, we need to distinguish the average Binder 
cumulant and its “global” counterpart ggi, depend¬ 
ing on when the disorder average is performed. They are 
defined as 


*?av — 


1 - 


3(m2)2_ 


dis 


5gi = 1 - 


[(m"^)]dis 

3[(m2)]L ■ 


( 6 ) 


The Binder cumulants g^^ and of the energy can be 
defined analogously. 

The correlation length is calculated via the second 
moment of the spin-spin correlation function G(r) = 
(1/L^) J2ij,a(^i‘Sf)S{r-rij),^^'^We again need to dis¬ 
tinguish average and “global” versions of this quantity, 
depending on when the disorder average is performed. 
They can be obtained efficiently from the Fourier trans¬ 
form G{q) of the correlation function: 

, ( 7 ) 

dis 

^ / [GjO) - G(gn.in)]dis \ 

Here, (jmin = 27r/L is the minimum wave number that 
fits into a system of linear size L. 

As was mentioned in Sec. lHI Al short Monte Carlo runs 
potentially introduce biases into observables for which a 
nonlinear operation is performed on the data before the 
disorder average. In our case, this includes C, y, ^av, and 


/ G(0)-G(ginin) y^' 

\ 9min^(9min) j 
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5av As is explained in Appendix these biases can be 
eliminated by using improved estimators. 

To judge the quality of the fits of our data to various 
mathematical models, we use the reduced weighted error 
sum x^. For fitting n data points {xi,yi) to a function 
f{x) containing q fit parameters, it is defined as 




1 f{xi)Y 

n — q ^ af 

I ^ 


(9) 


where af is the variance of yi. The fits are of good quality 

if X^ « 2 . 


IV. RESULTS 



A. Clean Ashkin-Teller model 

We first perform a number of simulations (using code 
A) of the clean model (no dilution, uniform interactions 
J = l,Ar = e) to test our algorithms and for later com¬ 
parison with the disordered case. 

For e = 0, the three-color Ashkin-Teller model is iden¬ 
tical to three independent Ising models. We simulate 
this model on lattices of 50^ to 800^ sites, averaging the 
data over 1000 samples for each size. The critical tem¬ 
perature is determined, as usual, from the crossing of the 
magnetic Binder cumulant g for different system sizes 
L. We find Tc = 2.26920(4) (the number in brackets 
indicates the error of the last digit) in agreement with 
the exact value 2/ln(l -|- 2 ^A) = 2.269185 .... Straight 
power-law fits (without subleading corrections) of mag¬ 
netization and susceptibility at Tc as functions of L yield 
the critical exponent estimates Plv = 0.1253(4) and 
xjv = 1.751(2). The correlation length exponent itself 
derives from the temperature derivative of g at critical¬ 
ity. We find v = 0.992(7). All hts are of good quality (re¬ 
duced of about 0.3 ... 1.2), and the estimates are in ex¬ 
cellent agreement with the exact values /3 = 1/8 ,7 = 7/4, 
and V = 1. 

For nonzero positive e, the phase transition in the clean 
Ashkin-Teller model is knowni^^— to be of first-order. 
We confirm this by simulations for e = 0.2, 0.3 and 0.5 
using systems of up to 560^ sites, averaged over 1000 
samples. Because of the supercritical slowing down as¬ 
sociated with hrst-order transitions we increase the equi¬ 
libration period to up to 500 Monte Carlo sweeps and 
the measurement period to 2000 sweeps. The Hrst-order 
character can be seen clearly in the double-peak structure 
of the probability distribution P{m) of the magnetization 
close to the transition temperature, as shown in Fig. [5] 
Notice that the minimum between the peaks becomes 
more pronounced with increasing system size, and the 
distance between the peaks remains roughly unchanged. 

We also perform exploratory simulations for e = 0.7 
and 1.0. The Wolff and Swendsen-Wang algorithms are 
invalid for these values because the effective interactions 
O can become negative. We therefore employ only 


FIG. 2. (Color online) Left: Magnetization distribution 
P(m) of the clean three-color Ashkin-Teller model for cou¬ 
pling e = 0.3 close to the transition temperature Tc ~ 3.178. 
The double-peak structure characteristic of a first-order tran¬ 
sition becomes more pronounced with increasing system size. 
(The curves for L = 400, 280, and 200 are shifted upwards by 
multiples of 0.5 for clarity.) Right: Wang-Landau density of 
states p{E) weighted by the Boltzmann factor (and normal¬ 
ized to its maximum) for the clean Ashkin-Teller model at 
e = 1. The system size is L = 48. 


Metropolis updates which requires long equilibration and 
measurement times and severely restricts the possible 
system sizes. Consequently, our simulations of up to 
200^ lattice sites (using 1000 samples, each with 5000 
equilibration sweeps and 10000 measurement sweeps) are 
less accurate than the simulations for e < 0.5. However, 
by comparing runs with “hot” and “cold” starts we can 
bracket the transition temperature with reasonable preci¬ 
sion. Analogous simulations of the clean system are also 
carried out using code B. 

The supercritical slowing down can be overcome by 
alternative sampling approaches 42r— To further check 
the correctness of the phase diagram, we therefore im¬ 
plement a code based on the Wang-Landau algorithm^ 
which is particularly suited to address first-order tran¬ 
sitions. This method performs a random walk in en¬ 
ergy space and provides direct access to the density of 
states p{E) (here, E = L?e refers to the extensive total 
energy). Specifically, the algorithm proceeds as follows: 
We initially set p{E) = 1. The energy histogram, which 
records the visit to each energy level E, is started at 
H{E) = 0. We then flip spins according to the proba¬ 
bility p{Ei —>■ Ej) = p{Ei)/p{Ej)) where Ei and 

Ej are the energies of the states before and after the 
flip, respectively. After every attempted spin flip, the 
density of states at the resulting energy is updated via 
p{E) -T fx p{E), and we record the visit by updating the 
histogram, H{E) —5> H{E) + 1. The modification factor / 
is initially set to / = exp(l). Once the histogram H{E) 
becomes reasonably flat, we reset H{E) = 0 and update 
the modification factor to a smaller value, / —>■ R- 
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FIG. 3. (Color online) Phase diagram of the clean three- 
color Ashkin-Teller model as function of temperature T and 
coupling strength e. The blue squares and the pink crosses 
mark the numerically determined transition points while the 
line is just a guide to the eye. The error bars of all our data are 
significantly smaller than the symbol size. For comparison, 
the figure also shows data extracted from the papers by Grest 
and Widoroi^ and Bellafard et ah— 

erating this procedure until / < exp( 10 “®) gives the den¬ 
sity of states p{E) with high precision. The right panel of 
Fig.[5]shows the resulting density of states, weighted with 
the Boltzmann factor, for the clean Ashkin-Teller model 
at e = 1. The double peak structure characteristic of two 
coexisting phases is clearly visible. The phase boundary 
can also be estimated from the peak of the specific heat 
curve. 

The phase diagram presented in Fig. [3] summarizes 
the results of our calculations for the clean Ashkin-Teller 
model. The figure also shows the critical temperatures 
reported in Refs. [H and[ll (extracted by redigitizing Fig. 
3 of Ref. [H and Fig. 1 of Ref. [T^. Within the errors of 
the redigitized data, our results agree well with Ref. [H 
but disagree with Ref. [l8l4^ 

B. Site-diluted Ising model 

After discussing the clean limit p = 0,e ^ 0, we now 
turn to the opposite limit p 7 ^ 0, e = 0. In this limit, our 
Hamiltonian is equivalent to three decoupled site-diluted 
Ising models. The critical behavior of the disordered two- 
dimensional Ising model is actually an interesting topic 
in itself because the clean correlation length exponent 
takes the value v = 1 which makes it marginal with re¬ 
spect to the Harris criterion^ dv > 2. In the literature, 
two main scenarios for the critical behavior have been 
put forward, the logarithmic correction scenario and the 
weak-universality scenario. 

The logarithmic correction (strong-universality) sce¬ 
nario arises from a perturbative renormalization-group 
approachi^^— It predicts that the asymptotic critical be¬ 


havior of the disordered Ising model is controlled by the 
clean Ising fixed point. Disorder, which is a marginally 
irrelevant operator, gives rise to universal logarithmic 
corrections to scaling. Specifically, one can derive the 
following finite-size scaling behaviori^— in the limit of 
large L. The specific heat at the critical temperature 
diverges as 

C ^ InlnL (10) 

with system size. Magnetization and magnetic suscepti¬ 
bility at Tc behave as 

M - [1 -k 0 (l/(lnL))] , ( 11 ) 

X~L^/ni + 0(l/(lnL))] , (12) 

with 7 /^ = 7/4 and fi/v = 1/8 as in the clean Ising 
model. Any quantity R of scale dimension zero (such as 
the Binder cumulants pav and g^i as well as the correla¬ 
tion length ratios ■Cav/D and ^gi/L) and its temperature 
derivative scale as 

R*+0{l/(lnL)) , (13) 

dR/dT - Li/^(lnL)-i/2 [1 + 0(l/(lnL))] (14) 

with ly = 1. This means, x, M, and R do not have 
multiplicative logarithmic corrections but dR/dT has a 
multiplicative (lnL)“^/^ correction. 

The weak-universality scenario was developed heuris- 
tically based on early numerical data42r— It states that 
the observables display simple power-law critical singu¬ 
larities. Their exponents vary continuously with disorder 
strength, but certain ratios stay constant at their clean 
values, for example ^/v and /d/v- The debate over the 
critical behavior of the two-dimensional disordered Ising 
model has persisted over many years, mainly because it is 
very hard to discriminate between logarithms and small 
powers on the basis of numerical data. Only in the last 
few years, the evidence seems to favor the logarithmic 
correction scenario (see, e.g., Refs.|43,lH,[5l, and Island 
references therein). 

The purpose of our simulations is twofold: On the one 
hand, the disordered Ising model is an important (lim¬ 
iting) reference case for our main topic, the disordered 
Ashkin-Teller model. On the other hand, we hope to 
make a contribution towards resolving the above contro¬ 
versy about the disordered Ising model itself. We there¬ 
fore perform a series of high-accuracy simulations for di¬ 
lution p = 0.3 and e = 0, using linear system sizes from 
L = 50 to 2240. The numbers of disorder realizations 
range from 4 x 10® for L = 50 to 5 x 10® for L = 2240). 

Figured] shows the Binder cumulant gg\ as a function 
of temperature. Analogous plots can be produced for gav, 
^gi/L and Cav/D. All these quantities display significant 
corrections to scaling manifest in the shift of the crossing 
temperature with increasing L. To extrapolate to infinite 
system size, we determine the crossings of the g^\ vs. T 
curves for sizes L/2 and L (and the corresponding cross¬ 
ings for (jfav, ^gi/L and ^av/F). Figure [3] presents the de¬ 
pendence of the crossing temperatures on the system 
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FIG. 4. (Color online) Binder cumulant g^\ vs. temperature 
T for p = 0.3 and e = 0 for different linear system sizes L. 
The statistical errors are much smaller than the symbol size. 
With increasing L, the crossing point shifts towards higher T, 
indicating significant corrections to scaling. 



FIG. 5. (Color online) Crossing temperatures vs. inverse 
system size 1/L for p = 0.3 and e = 0. Tr is the temperature 
where the curves of pav, Pgi, Cav/T and ^gi/L versus T cross 
for system sizes Lj2 and L. The solid lines are fits to Tc = 
Tc + aL~^. The error bars of are about the size of the 
symbols at the right side of the plot and become much smaller 
towards the left. 


size. The critical temperature Tc can be extracted by 
extrapolating the crossing temperatures to infinite sys¬ 
tem size. Fits to = Tc + aL~^ yield Tc = 1.07201(3) 
which agrees reasonably well with the result of Ref. 1^, 
Tc = 1.07194(6), obtained from systems with up to 256^ 
sites. (Fits of T^, vs. 1/T to quadratic polynomials give 
comparable results.) 


1. Logarithmic-correction scenario 

To study the critical behavior, we analyze the system- 
size dependence at the critical temperature of magne¬ 
tization M, susceptibility y, specific heat C and the 
slope filn(^gi/T)/fiT of the normalized correlation length 
curves. Straight power-law fits (without corrections to 
scaling) of the data at T = 1.07200 to M ^ 

X C - T“/'^ and d\n{^gi/L)/dT - L^/'' give 

the estimates ^jv = 0.1217(1), y/jz = 1.8046(2), ajv = 
0.0516(1), and v = 1.107(3). These values do not agree 
with the clean Ising exponents, /3/iz = 1/8, 7 / 1 / = 7/4, 
ajv = 0, and v = 1. However, the quality of all fits 
is extremely poor, with reduced values of about 50, 
1300, 4600, and 7, respectively. This indicates that the 
data significantly deviate from pure power laws. 

To understand the nature of the deviations, we divide 
out the clean Ising power laws and plot the resulting data 
in Fig.ini MT^/® and shown in the upper panel, 

clearly increase much more slowly than power laws with 
L. As suggested in eqs. (ED and ED, their behaviors 
can be analyzed as logarithmic corrections to scaling and 
fitted to the form a[l -I- 6 /ln(cT)]. The fits are of good 
quality (reduced y^ of 1.1 and 0.5 for the magnetization 
and susceptibility, respectively). The lower panel of Fig. 
m shows specific heat C vs. system size at the critical 
temperature. The data can be fitted well by the double- 
logarithmic form aln[61n(cT)] suggested by eq. (fTUl) . giv¬ 
ing a reduced y^ of about 1.2. For comparison, we also 
consider a simple logarithmic form C = a ln( 6 T). A semi¬ 
log plot of C vs. ln(T) (not shown) shows strong devia¬ 
tions from a straight line. Correspondingly, the simple 
logarithmic fit is of very poor quality, with a reduced y^ 
of about 3000; and it does not improve much if the fit 
range is restricted. 

The lower panel of Fig. [ 6 ] also shows the slopes 
dln(^av/T)/dT and dln(^gi/T)/dT of the normalized 
correlation length curves at the critical tempera¬ 
ture. We again divide out the clean Ising power law 
d ln(^/T)/ dT ^ T to make the deviations from power-law 
behavior clearly visible. The resulting data can be fitted 
well by the logarithmic form a[ln( 6 T)]“^/^ suggested by 
eq. (|TT)) (reduced y^ of about 0.3 for both data sets). 
Including extra additive corrections to scaling does not 
improve the fits. The slopes of the magnetic Binder cu- 
mulants behave analogously. 

In addition to the quantities shown in Fig.[6l we study 
the system-size dependence of the dimensionless ratio 
f/T at criticality. Within the logarithmic correction sce¬ 
nario, this ratio is expected to approach the universal 
value (^/T)* of the clean Ising model which is known with 
high precision (see, e.g.. Ref. [13). For a square lattice 
with periodic boundary conditions (torus topology), it 
reads (^/T)* = 0.9050488292(4). According to eq. (fT5]) . 
the approach to this value is logarithmically slow. We 
therefore plot ^gi/T and ^av/T at the critical tempera¬ 
ture as functions of l/ln(T) in Fig. [71 Both ratios can 
be well fitted to the form (^/T) = (^/T)* -f a/ln(6T), 
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FIG. 6. (Color online) System-size dependence of observables 
at the critical temperature. Top: Log-log plots of and 

^^-7/4 L at T = 1.07200 for p = 0.3 and e = 0. The solid 
lines are hts to a[l -|- 6/ln(cL)] as suggested by eqs. (Illll and 
m- Bottom: Log-log plots of the specific heat C as well as 
the slopes L“^dln(^av/L)/dT and L)/dT vs. L at 

T — 1.07200 for p = 0.3 and e = 0. The solid lines represent 
fits to aln[61n(cl/)] for the specihc heat and to a[ln(6I/)]“^/^ 
for the slopes. The dashed line shows a power-law fit using 
V = 1.130 as implied by the hyperscaling relation 2 — a = 2v 
in the weak-universality scenario. 


as suggested by eq. (USD, with {^/L)* fixed at the clean 
Ising value. The fits are of excellent quality (reduced 
of 0.7 and 0.6, respectively) if the fit range is re¬ 
stricted to system sizes L > 70. We attribute the small 
deviations for the smallest L to subleading terms^ of 
the form lnln(6L)/ln^(&L) in eq. (ITSl) that are not in¬ 
cluded in the fit. Analyzing the system size dependence 
of the Binder cumulant g^-v at criticality gives the same 
result: (;av approaches the universal clean Ising valued, 
g* = 0.610692(2) following eq. (IT^ . The behavior of g^x 
is more complex. With increasing L, it first decreases 
below g* before turning around and approaching g* from 
below. A quantitative analysis therefore requires even 
larger systems than ours to properly fit the subleading 
terms in eq. 0. 


FIG. 7. (Color online) Dimensionless ratios 5gi/L and ^av/L 
at criticality vs. l/\n{L) for p = 0.3, e = 0 (T = 1.07200), 
p = 0.3, e = 0.5 (T = 1.93471), and for p = 0.1, e = 0.1 (T = 
2.18769). The error bars are significantly smaller than the 
symbol size. The lines are fits to (^/L) = (C/T)* + a/ln(6L) 
with i^/L)* fixed at the clean Ising value 0.9050488292. Note 
that the two ^gi/L curves for p = 0.3, e = 0 and p — 0.3, e = 
0.5 are almost on top of each other. 


Finally, we also study the product order parameter Mp. 
As the different colors are completely independent for 
e = 0, Mp must scale as Indeed, the system size 
dependence of our data at criticality (not shown) can be 
fitted very well by Mp = -I- b/\n{cL)]. 

In sum, our high-accuracy data almost perfectly agree 
with the renormalization-group predictions that lead to 
the logarithmic correction scenario outlined in eqs. (nni) 
to (fTH) over the entire range of system sizes studied {L = 
50 to 2240). 


2. Weak-universality scenario 

Can these data also be understood within the heuris¬ 
tic weak-universality scenario? Figure [5] shows that the 
data deviate significantly from pure power laws over en¬ 
tire system size range. The weak-universality scenario 
can thus only work, if at all, if corrections to scaling are 
included (in addition to potential changes in the critical 
exponents). 

The system size dependencies of M and x at Tc shown 
in Fig. [5] can be fitted with the clean Ising exponents 
Pjv = 1/8 and y/zz = 1, provided that corrections 
to scaling of the type M = -|- bL~^) and 

X = bL~^) are included. These fits are of 

lower, but still acceptable, quality (reduced x^ of about 
1.9 and 2.5, respectively) than the fits with logarithmic 
corrections given in eqs. (HU and (HU- Four-parameter 
fits to the same functional forms but with floating crit¬ 
ical exponents give Pjv = 0.123(10) and y/zz = 1.76(2) 
where the errors mostly stem from the sensitivity of the 
fits towards changes of the fit interval. We conclude that 
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(ijv and ^jv agree with the clean Ising values. As these 
exponent ratios are expected to take the clean values in 
both scenarios, this does not allow us to discriminate be¬ 
tween the scenarios. We therefore turn to the exponents 
ajv and v which are expected to be nonuniversal. 

In the weak-universality scenario, the slow increase of 
the specific heat C with L (as shown in the lower panel of 
Fig. [H) is interpreted as power-law behavior of the type 
C = Coo + with a negative exponent a. A fit 

to this form yields ajv = —0.230(3), however, it is of 
much lower quality (reduced of about 26) than the 
double-logarithmic fit employed above. Moreover, the fit 
is very unstable. If we extract an effective exponent by 
restricting the fit to the interval (Amim 4Li„in)i its value 
increases monotonically from —0.278 for Amin = 50 to 
—0.106 for Amin = 560, with no sign of saturation. A 
controlled extrapolation to infinite Amin is difficult; but 
it appears to be compatible with ajv = 0. In fact, a 
(perhaps overambitious) five-parameter fit that includes 
subleading corrections to scaling, C = Coo + aA“/'^(l -|- 
bL~‘^), yields a very small ajv ^ 0.03 albeit with a large 
error of about 0 . 2 . 

The temperature derivatives dln(^av/A)/dr and 
(iln(^gi/A)/(iT of the normalized correlation length 
curves at the critical temperature can be fitted with the 
clean Ising exponent v = \ \i corrections to scaling of 
the type aA^/‘"(l -|- bL~^) are included. The quality of 
these fits (reduced of 0.5 and 0.2) is comparable to 
that of the logarithmic fits above (note, however, that 
the logarithmic fits contained only two free parameters). 
Four-parameter fits to the same functional form but with 
floating v give v = 1.07(4) and 1.06(6) for the average 
and global correlation length data, respectively. 

If we ignore the strong size-dependence of the effec¬ 
tive specific heat exponent and take the value ajv = 
—0.230(3) resulting from the global fit, the hyperscaling 
relation 2 — a = 2v yields a correlation length exponent 
of zz = 1.130(3). As the lower panel of Fig. | 6 ] shows, 
the d\w{^lL)ldT data are clearly incompatible with this 
value, even if corrections to scaling are includedi^ 

Finally, we point out that the subleading exponent uj 
appearing in all the power-law fits is not very robust. Its 
values seem to cluster around 0.35 but they vary between 
about 0.2 and 0.7 upon changing the fit intervals and 
between quantities. 

We conclude that the weak-universality scenario is not 
compatible with our numerical results: Simple power-law 
singularities do not describe the data at all. If corrections 
to scaling are included, we do not find evidence for the 
asymptotic exponents to be different from the clean Ising 
ones. 


3. Other dilutions 

We have performed analogous simulations for dilution 
p = 0.2, using systems with 50^ to 1120^ sites. The 
data are averaged over 10® to 10® disorder configurations. 



FIG. 8. (Color online) Magnetization distribution P{rn) of 
the three-color Ashkin-Teller model with dilution p = 0.3 and 
coupling e = 0.5 close to the transition temperature Tc « 
1.94. The distribution features a single peak characteristic of 
a continuous transition. (The curves for A = 400, 280, and 
200 are shifted upwards by multiples of 1.0 for clarity.) 


By extrapolation the crossing temperatures of the Binder 
cumulant and the normalized correlation length as above, 
we find a critical temperature of = 1.50709(5). The 
system-size dependence of observables at Tc looks almost 
identical to that shown in Fig. [6]for p = 0.3: The data 
feature pronounced deviations from power-law behavior 
that can be fitted very well by the logarithmic-correction 
scenario, eqs. (|T(I1) to (|TH) . over the entire range of system 
sizes studied (the reduced range between 0.4 and 1.1). 


C. Site-diluted Ashkin-Teller model 

After having discussed the limiting cases, we now turn 
to the full problem, the site-diluted Ashkin-Teller model. 
We first consider a system with dilution p = 0.3 and four- 
spin coupling e = 0.5 because we expect deviations from 
the clean Ising critical behavior, if any, to be more easily 
visible if p and e are large. 

According to the Aizenman-Wehr theorem^ the first- 
order phase transition of the clean system should be de¬ 
stroyed by dilution. We confirm this by calculating the 
magnetization distribution P{m) close to the transition 
temperature for systems of 200^ to 560^ sites (1000 dis¬ 
order configurations each). It is shown in Fig. |8l The 
distribution features a single broad peak characteristic 
of a continuous transition, in contrast to the double-peak 
structure of the clean case (Fig. [ 5 )). 

We then perform a series of high-accuracy simulations 
for linear system sizes from A = 50 to 1600. The num¬ 
bers of disorder configurations range from 3 x 10® for 
A = 50 to 5 X 10® for A = 1600. To find the critical 
point, we study the Binder cumulants Pav(A) and 5 gi(A) 
as well as the normalized correlation lengths ^av(A)/A 
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FIG. 9. (Color online) Normalized correlation length £,^\/L 
vs. temperature T for p = 0.3 and e = 0.5 for different linear 
system sizes L. The statistical errors are signihcantly smaller 
than the symbol size. With increasing L, the crossing point 
shifts towards lower T, indicating significant corrections to 
scaling. 



1/L 

FIG. 10. (Color online) Crossing temperatures Tx vs. inverse 
system size 1/L for p = 0.3 and e = 0.5. Tx is the temperature 
where the curves of gav, ( 7 gi, Cav/L and ^gi/L versus T cross 
for system sizes L/2 and L. The solid lines are hts to Tx = 
Tc + aL~^. The error bars of Tx are about the size of the 
symbols at the right side of the plot and become much smaller 
towards the left. 


and ^gi(T)/L. The resulting data look qualitatively very 
similar to those of the diluted Ising model. As an exam¬ 
ple, we present ^gi/L in Fig.[3 The critical temperature 
can be estimated by extrapolating to infinite system size 
the temperatures where the ga.v{T) curves (as well as the 
ggi{T), ^av(F)/F and ^g\{T)/L curves) for sizes L/2 and 
L cross. Figure [To] shows the system-size dependence of 
the crossing temperatures T^- Fits to + aL~^ 

yield the estimate = 1.93472(5). (Fits of vs. 1/L 
to quadratic polynomials give comparable results.) 


1. Critical behavior: Ising with logarithmic corrections 

To analyze the critical behavior, we now study the 
system-size dependence at criticality of magnetization 
M, susceptibility y, specific heat C, and the slope 
dln(^gi/L)/dr of the normalized correlation length. Sim¬ 
ple power-law fits over the entire system size range (L = 
50 to 1600) of the data at L = 1.93471 to M 
X L'^/'^, C - L“/'^ and dln{^gi/L)/dT - L^/'' give 
the estimates jdjv = 0.1238(1), y/iz = 1.7948(3), ajv = 
0.0735(1), and v = 1.075(2). These values do not agree 
with the clean Ising exponents, but the quality of the 
fits is again very poor. The reduced values are about 
21, 180, 4300, and 7, respectively, indicating systematic 
deviations from pure power-law behavior. 

To investigate these deviations in detail, we proceed 
analogously to the diluted Ising model in Sec. IIVBI i.e., 
we divide out the clean Ising critical behavior and present 
the resulting data in Fig.[TTJ The figure shows that none 
of the plotted quantities follow simple power laws; instead 
they vary more slowly with L over the entire system size 
range. 

Motivated by Cardy’s renormalization group^^, we 
therefore attempt to fit the data with the clean Ising ex¬ 
ponents and logarithmic corrections analogous to those 
of the diluted Ising model, eqs. m to m- The magne¬ 
tization can be fitted well to the form M = aL“^/®[l -|- 
hj In(cL)] over the entire size range L = 50 to 1600. The 
reduced y^ is about 0.9. Fitting the susceptibility to 
y = aL^/^[l -I- 5/ln(cL)] over the entire size range leads 
to an unsatisfactory reduced y^ ss 7. However, the fit be¬ 
comes of good quality (reduced y^ ~ 1.8) if we drop the 
two smallest system sizes, restricting the fit to the range 
L = 100 to 1600. We attribute this to the crossover from 
the strong first-order transition in the clean case to our 
critical point. (This crossover will be studied in detail in 
Sec, live 31 ) 

We also analyze the product order parameter Mp. 
Within Cardy’s theory^ the critical renormalization 
group fixed point is at e = 0. This means different colors 
decouple at criticality. Thus, Mp should scale as M^. In 
agreement with this expectation, the system size depen¬ 
dence of our data at criticality (shown in the inset) can 
be fitted by Mp = aL“^/‘^[l-|-6/ln(cL)], giving a reduced 
y^ of about 1. 

The specific heat, shown in the lower panel of Fig. |Tl] 
can be fitted well by the double-logarithmic form C = 
aln[61n(cL)] over the entire size range, giving a reduced 
y^ of a about 1.2. Finally, the slopes dln(^av/L)/dr and 
dln(^gi/L)/dr of the normalized correlation lengths at 
criticality can be fitted by the form aL[ln(6L)]“^/^ over 
the entire size range (reduced y^ of about 0.5 and 0.2, 
respectively). 

Finally, we investigate the system size dependence of 
the dimensionless ratio ^/L at criticality. Fig. |7| presents 
^gi/L and ^av/L as functions of 1/ ln(L). Both ratios can 
be well fitted to the form (^/L) = {^/L)* + a/\n{bL) 
with (^/L)* fixed at the clean Ising value, as suggested 
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FIG. 11. (Color online) System-size dependence of observ¬ 
ables at the critical temperature. Top: Log-log plots of 
and vs. L at T = 1.93471 for p = 0.3 and 

e = 0.5. Inset: Log-log plot of vs. L. All solid 

lines are fits to a[l -|- b/ln{cL)]. Bottom: Log-log plots of 
the specific heat C as well as the slopes L~^/L)/dT 
and Z/“^dln(^gi/L)/dT vs. L a.t T = 1.93471 for p = 0.3 and 
e = 0. The solid lines represent fits to aln[61n(cL)] for the 
specific heat and to a[ln(feL)]“^/^ for the slopes. 


by eq. (USD. The fits are of good quality (reduced of 
1.0 and 1.6, respectively) if the fit range is restricted to 
system sizes L > 100. The deviations for the smaller L 
likely stem from the crossover between the clean first- 
order transition and our critical point as well as from 
subleading terms of the form In ln(6L)/ In^ (bL) in eq. (fT^ 
that are not included in the fit. 


We conclude that all our data can be described nearly 
perfectly in terms of the clean Ising critical behavior 
with logarithmic corrections to scaling, as predicted by 
Cardy’s renormalization group:^ 


2. Power law behavior? 

Even though our analysis does not show any disagree¬ 
ments between the Monte Carlo data and renormaliza¬ 
tion group predictions, we still test whether the data are 
compatible with nonuniversal power-law critical behav¬ 
ior as suggested in Refs. 0 and 0 Since the quantities 
shown in Fig. [TT] do not follow simple power laws, it is 
clear that corrections to scaling need to be included in 
addition to possible deviations of the exponents from the 
clean Ising values. 

Magnetization, susceptibility and the slopes of the 
normalized correlation lengths can be fitted to M = 
aL-/3/i'(l + bL-‘^), X = aLT'/‘'(l -b and 

d\n{£_/L)/dT = + bL~‘^). If the exponents/3/i^ 

and y/zz and v are fixed at the clean Ising values 1/8, 
7/4 and 1, respectively, the fits are of good quality, with 
reduced just slightly higher than the logarithmic fits 
above. (The system size range for the susceptibility fit 
needs to be restricted to L > 100 to achieve an accept¬ 
able quality.) Four-parameter fits over the entire system 
size range to the same functional forms, but with float¬ 
ing (djv, y/z/, and v give the values /3/zz = 0.125(1) and 
y/zz = 1.78(1) and v = 1.04(6) where the errors mostly 
stem from the sensitivity of the fits towards removing 
points from the ends of the system size range. Note that 
Pjv and y/zz do not quite fulfill the hyperscaling relation 
2/3/zz -b y/zz = 2 , suggesting that these values are not the 
true asymptotic exponents. 

It is worth pointing out that the effective inverse cor¬ 
relation length exponent 1 /zzeff, obtained by fitting the 
correlation length slopes over a finite system size range, 
is always smaller than unity. This can be seen from the 
downward slope of L~^d\n{^/L)/dT vs. L in the lower 
panel of Fig. [TT] it also directly follows from (fTTl) . The 
effective correlation length exponent thus fulfills zzgff > 1 . 
In contrast to Refs. 0 an d 0 we see no indications of the 
inequality dv > 2 due to Chayes et alj^ being violated 
even by the effective exponent. 

The most interesting quantity is the specific heat. A 
correlation length exponent zz > 1 implies, via the hy¬ 
perscaling relation 2 — a = 2 zz, that the specific heat 
exponent a < 0. We therefore attempt to fit the data to 
the form C = Coo + clL°‘^'^ . A ht of all system sizes yields 
aju = —0.170(4), but the reduced ~ 12 is unaccept¬ 
ably large. Moreover, the ht is unstable. If we extract 
an effective exponent (Q;/zz)egf by restricting the ht in¬ 
terval to (Amin, 4L„iin)) its value varies between —0.204 
and —0.131 for Amin between 50 and 400. Extrapolating 
these values to inhnite system size does not give a deh- 
nite answer. Depending on the mathematical model used 
for the extrapolation, we hnd values between — 0.12 and 
0 . 

We again point out that the subleading exponent w 
appearing in all the power-law hts is not very robust. Its 
values vary between about 0.2 and 1.1 upon changing the 
ht intervals and between quantities. 

We conclude that the description of our data in terms 
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FIG. 12. (Color online) Semi-log plot of specific heat C vs. 
system size L at criticality for dilution p = 0.1 and several 
couplings e. The error bars are much smaller than the symbol 
size. The solid lines are fits to C = aln[61n(cL)]. 


of power-law singularities does not work nearly as well as 
the logarithmic correction scenario of Sec. IIV C II Sim¬ 
ple power laws do not describe the data. If we insist on 
fitting the data to power laws with corrections to scal¬ 
ing included, there is no compelling evidence for the true 
asymptotic exponents to differ from the clean Ising val¬ 
ues. 


3. Universality and crossover between the elean and dirty 
phase transitions 

We perform analogous simulations for several addi¬ 
tional values of dilution p and coupling strength e in 
order to test whether the asymptotic critical behavior 
is universal. In addition, we wish to explore the interest¬ 
ing crossover from the first-order transition of the clean 
Ashkin-Teller model to the critical point of the disordered 
system. As discussed at the end of Sec. HTJ it should be 
particularly pronounced when the first-order transition 
of the clean system is strong and the disorder is weak. 

To analyze the crossover, we therefore perform a se¬ 
ries of simulations for the weaker dilution p = 0.1. The 
coupling e takes values 0.1, 0.2, 0.3, and 0.5. (Increasing 
e increases the strength of the first-order transition in 
the corresponding clean system.) These simulations use 
sizes between L = 25 and 1120 with 10^ to 10® disorder 
realizations each. The data analysis follows the steps out¬ 
lined above, and the resulting critical temperatures are 
listed in the legend of Fig. [T^l Interestingly, the crossing 
temperature of the Binder cumulants and the correla¬ 
tion length ratios shifts much less with system size than 
for p — 0.3, indicating weaker disorder-induced correc¬ 
tions to scaling. (For the crossings between the L = 50 
and L = 100 curves, {T^ — Tc)/Tc is roughly one order of 
magnitude smaller for p = 0.1 than for p = 0.3.) 



FIG. 13. (Color online) Semi-log plot of ys. L at 

criticality for dilution p — 0.1 and several couplings e. The 
data for e = 0.2, 0.3 and 0.5 are shifted upwards by 0.003, 
0.005, and 0.008 for clarity. The error bars are much smaller 
than the symbol size. The solid lines are fits to = 

a[l + 6/ln(cL)], the fit ranges are indicated in the graph. 


Figure [12] displays a semi-logarithmic plot of the spe¬ 
cific heat C at criticality vs. system size L. For all e, C 
curves downward, indicating that it increases more slowly 
than logarithmic with L. The figure also shows fits to 
the double-logarithmic form C=aln[61n(cL)] suggested 
by m- While the fits look nearly perfect to the eye, a 
analysis reveals the effects of the crossover from clean 
to dirty behavior: For e = 0.1 and 0.2, fits over the entire 
size-range L = 25 to 1120 are of good quality (reduced 
« 0.4 and 0.7, respectively). The quality decreases for 
e = 0.3 (reduced « 2.3) and e = 0.5 (reduced ~ 9). 
Good quality fits (reduced y^ < 2) can be restored by 
restricting the fit range to L > 35 for e = 0.3 and to 
L > 50 for e = 0.5. 

More pronounced signatures of the crossover from 
clean to dirty behavior can be found in the corrections 
to the leading power-law size dependencies of magneti¬ 
zation and susceptibility (probably because these correc¬ 
tions are weak — they only change the observables by a 
few percent over the entire system-size range). Figure [T51 
presents y/L^/^ vs. L at criticality for dilution p = 0.1 
and several e. The data for e = 0.1 behave analogously 
to those observed earlier in Fig. [TT] forp = 0.3, e = 0.5. 
They can be fitted well (reduced y^ ~ 0.5) with the log¬ 
arithmic form y/L^/'^ = a[l -|- b/ In(cL)]. The same holds 
for the e = 0.2 data which give a reduced y^ ~ 0.7 (Note, 
however, that the curvature of the e = 0.2 data is very 
weak.) For larger e, the behavior changes. first 

decreases with increasing L before turning around and 
starting to increase slowly. The increase can be fitted 
to a[l -I- b/ In(cL)] for L > 200 (e = 0.3) and L > 400 
(e = 0.5). 

The corrections to the leading power-law behavior of 
the magnetization behave in the same fashion as those of 
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the susceptibility. We thus conclude that the asymptotic 
critical behavior is compatible with the logarithmic cor¬ 
rection scenario for all shown e values. We emphasize, 
however, that large system sizes are necessary to reach 
this asymptotic behavior if the first-order transition of 
the corresponding clean system is strong (e = 0.3 and 
0.5) and the disorder is weak. The universality of the 
critical behavior is also confirmed by the analysis of the 
ratio L at criticality for p = 0.1, e = 0.1, shown in Fig. 

m 

In addition to the parameter sets already discussed, 
we also perform simulations for p = 0.05, e = 0.05 as well 
as p = 0.3, e = 0.3 (system sizes between L = 50 and 
1120 with up to 10® disorder realizations). In both cases, 
the critical behavior can be fitted well with clean Ising 
critical behavior and logarithmic corrections to scaling 
over the entire system size range. 


D. Random-bond Ashkin-Teller model 

The random-bond Ashkin-Teller model is expected to 
be in the same universality class as the site-diluted model 
because both types of randomness are implementations of 
random-Tc disorder. However, in view of the unexpected 
results of Refs. in and n we also perform a number of 
simulations for the random-bond case. 

We employ code A to simulate systems using the bi¬ 
nary bond distribution ([2|) with Jh = 2, J; = 0.5 and 
concentration c = 0.5. The four-spin interactions Kij are 
slaved to the Ising interactions Jij via Kij = eJij with 
uniform e. We perform a series of runs for e = 0,0.1,0.2, 
and 0.5. The linear system sizes are between L = 35 and 
1120 ; and the numbers of disorder realizations for each 
parameter set range from 10® for L = 1120 to 10® for 
L = 35. 

The data analysis follows the steps outlined in Sec. 
IIV Cl The critical temperatures found by extrapolating 
the crossing temperatures of the Binder cumulants pav 
and pgi as well as the correlation lengths ratios ^av/L and 
^gi/L to infinite system size are shown in the legend of 
Fig. [TT] We note that even though the bond random¬ 
ness looks substantial {Jh/Ji = 4), the disorder-induced 
corrections to scaling turn out to be rather weak: The 
shifts of the crossing temperatures with system size 
are even smaller than those for dilution p = 0.1. Because 
of the weaker corrections to scaling, effective exponents 
extracted by simple power-law fits over the entire sys¬ 
tem size range are already very close to the expected 
clean Ising values. For example, for e = 0.2, we find 
15jv = 0.1258(1), 7 / 1 / = 1.7527(5), and 1 / = 1.031(2). 

We now analyze whether the asymptotic critical be¬ 
havior can be described by logarithmic corrections to the 
clean Ising power laws, as was the case for site dilution. 
Figure [M] displays a semi-logarithmic plot of the spe¬ 
cific heat C at criticality vs. system size L. For all e, C 
curves downward, indicating that it increases more slowly 
than logarithmic with L. The figure also shows fits to 



FIG. 14. (Color online) Semi-log plot of specific heat C vs. 
system size L at criticality for the random bond-Ashkin-Teller 
model with Jh = 2, Ji = 0.5 and c = 0.5. The error bars are 
much smaller than the symbol size. The solid lines are fits to 
C = aln[61n(cZ/)]. 



FIG. 15. (Color online) Semi-log plot of ys. L at 

criticality for the random bond-Ashkin-Teller model with 
Jh = 2, J( = 0.5 and c = 0.5. The data for t = 0.1, 0.2 and 
0.5 are shifted upwards by 0.004, 0.007, and 0.013 for clarity. 
The error bars are much smaller than the symbol size. The 
solid lines are fits to = a[l -|- 6/ln(cI/)], the fit ranges 

are indicated in the graph. 


the double-logarithmic form C=aln[61n(cL)] suggested 
by cni). All fits are of high quality (reduced < 2), 
for e = 0.5 this requires us to restrict the fit range to 
L > 100. 

The analysis of susceptibility and magnetization at 
criticality again reveal signatures of the crossover from 
clean to dirty behavior. Figure [TSl presents vs. 

L at criticality for different e. The data for e = 0 and 
0.1 can be fitted to the logarithmic form = a[l -I- 

6 /ln(cL)] for sizes L > 100 (with reduced y^ < 2). For 
larger e, first shows a pronounced decrease with 
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FIG. 16. (Color online) Log-log plot of the slopes 
dln(^gi/L)/cir vs. L at criticality for the random bond- 
Ashkin-Teller model with Jh = 2, Ji = 0.5 and c = 0.5. 
The data for e = 0.1, 0.2 and 0.5 are multiplied by factors 
2, 4, and 8 for clarity. The error bars are much smaller than 
the symbol size. The solid lines are hts to d\n(^^\/L)/dT = 
aL[ln(6L)]-i/^ 


increasing L before turning around and starting to in¬ 
crease slowly. The increase can be fitted to a[l+h/ In(cL)] 
for L > 200 (e = 0.2) and L > 280 (e = 0.5). It must 
be noted however, that the susceptibility corrections to 
scaling in these data sets are so weak (in agreement with 
the small shifts of mentioned above) that we can¬ 
not unequivocally confirm their functional form. Their 
large-L behavior is certainly compatible with the pre¬ 
dicted a[l -b 6 /ln(cL)] form but other functions would 
work as well. The corrections to the clean Ising power 
laws for the magnetization behave analogously to those 
of the susceptibility. 

Finally, Fig. fTHl shows the slopes /L)/dT of the 

correlation length ratio vs. L at criticality. All curves can 
be fitted with high quality (reduced ~ 1 ) by the form 
aL[ln( 6 L)]“^/^ suggested by eq. (fTH) . (Because the devia¬ 
tions from the clean Ising power laws are again weak, the 
fits cannot unambiguously discriminate between different 
functional forms: Simple power laws work as well, giving 
exponents v in the range of 1.03 to 1.05.) We conclude 
that the asymptotic critical behavior of the random-bond 
Ashkin-Teller model is fully compatible with Cardy’s pre¬ 
dictions, i.e., clean Ising exponents with logarithmic cor¬ 
rections. 

In all of the above (code A) simulations, the four-spin 
interactions Kij are slaved to the Ising interactions Jij 
via Kij = eJij with uniform e. In addition, we study 
random-bond Ashkin-Teller models with uniform Kij = 
K employing code B. As summarized in Sec. IIII Cl we 
consider Jh = 6/5 and J; = 4/5 with equal probability 
c = 0.5, implying an average interaction oi J = {Jh + 
Ji)/2 = 1. Note that this disorder is much weaker than 
the disorder considered in the code-A simulations above. 



FIG. 17. (Color online) Log-log plot of the magnetization M, 
susceptibility x, polarization Mp and polarization susceptibil¬ 
ity Xp vs. L at Tc for the random bond-Ashkin-Teller model 
with Jh = 6/5, Ji = 4/5 and uniform K — 0.1. All error bars 
are much smaller than the symbol sizes. The solid lines are 
power-law fits. 


The uniform, non-random four-spin interaction is given 
hy K = eJ with e = 0.1. We simulate systems having 
linear sizes from L = 24 to 1600. The number of disorder 
realizations ranges from 10^ for L = 1600 to 10® for L = 
24. 

The analysis of the data generated by code B proceeds 
as in Sec. IIV Cl The critical temperature is found by ex¬ 
trapolating the crossing temperature of the Binder 
cumulants gav to the infinite system size limit. This 
yields a critical temperature Tc = 2.55625(1). Simple 
power-law fits, shown in Fig. [T71 of observables at to 
M ~ X - Mp - and Xp ^ 

give PIv = 0.125(1), yjv = 1.749(3), (ip/v = 0.230(1) 
and Xp/v = 1.53(1). The fits are of good quality (once 
again reduced < 2) if we restrict them to system 
sizes L > 96. The exponents of magnetization and 
susceptibility have already locked onto the clean Ising 
values fdjv = 1/8 and yjv = 7/4 within their error 
bars. The exponents related to Mp and Xp do not quite 
agree with the expected values jdp/u = ^jdjv = 1/4 and 
Xp/v = 2 — 2j5p/v = 3/2, but they are close. This is in 
tune with the results obtained for random-bond Ashkin- 
Teller model simulated via code A. 

Can the deviations of Mp and Xp from the expected 
behavior be explained by logarithmic corrections? To 
answer this question, we again divide out the expected 
power laws and present the resulting data in Fig. [THl 
The product order parameter, Mp and the associated 
susceptibility Xp can be fitted quite well with Mp = 
q,7^-i/4(i-I- 5/in(cL)] and Xp = aL®/^[l-b 5/In(cL)]. For 
sizes L > 96, the reduced x^ are 1.26 and 0.713 for Mp, 
and Xp) respectively. 

Corrections to the clean Ising behavior of magnetiza¬ 
tion and susceptibility are very weak (in agreement with 
the fact that the exponents of simple power-law fits al- 
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FIG. 18. (Color online) Semi-log plot of and 

vs. L at Tc for the random bond-Ashkin-Teller model with 
Jh = 6/5, Ji = 4/5 and uniform K = 0.1. The solid lines are 
fits to a[l -I- b/ln{cL)]. Insert: Semi-log plot of the specific 
heat C vs. L, the solid line is the fit to aln{bL). 

ready coincide with the clean Ising ones). If we include 
the smaller system sizes in the fits, these weak correc¬ 
tions cannot be fitted satisfactorily with the universal 
logarithms (El) and ED- Similarly, the specific heat does 
not follow the double-logarithmic form, C = oln[61n(cL)] 
(see inset of Fig.[TH]). Instead, the data for large system 
sizes L > 768 are best described by the single logarithmic 
form, C = a In (bL) expected at the clean Ising critical 
point. 

How can we explain these observations? In the present 
system, both the bare disorder strength and the coupling 
e are rather weak. The renormalization group flow (Fig. 
[T|) therefore does not travel too far from the origin, ex¬ 
plaining that our effective exponents are very close to 
the clean Ising ones. Note, however, that the renormal¬ 
ization group “time” (flow parameter) and, correspond¬ 
ingly, the system size range needed to go around the loop 
in Fig. [1] do not vary much with the size of the loop. 
As the bare disorder is weak, the system thus does not 
reach the falling, asymptotic part of the loop even for 
L = 1600. This may also explain why the effective cor¬ 
relation length exponent = 0.93(2) that we extract 
from the slope of the Binder cumulant vs. temperature 
curves in this weakly disordered system does not fulfill 
the Chayes’ inequality^ dv > 2. Because of the expo¬ 
nential dependence of the breakup length on the dis¬ 
order strength, confirming this picture numerically would 
likely require enormous system sizes. 

V. CONCLUSIONS 

In summary, we have performed high-accuracy Monte 
Carlo simulations of the disordered three-color Ashkin- 
Teller model in two dimensions using systems with up 
to 1600^ lattice sites. We have investigated two types 
of disorder, random site dilution and random interac¬ 


tions (bond randomness). Our results show that the first- 
order phase transition of the clean Ashkin-Teller model 
is destroyed by the randomness, in agreement with the 
Aizenman-Wehr theorem»ii^ 

We have carefully analyzed the critical behavior of the 
h emerging continuous phase transition and found strong 
evidence that the asymptotic critical behavior is univer¬ 
sal and agrees with the predictions of Cardy’s renormal¬ 
ization group theoryiii This means, the critical expo¬ 
nents coincide with those of the clean two-dimensional 
Ising model, but with additional logarithmic corrections 
to scaling analogous to those found in the disordered two- 
dimensional Ising model. For example, the specific heat 
takes the characteristic double-logarithmic form ED- 
What could be the reason for the differences between 
our results and the unusual behavior (nonuniversal criti¬ 
cal exponents and violations of the inequality dv >2 due 
to Chayes et ali^) reported in Refs. [3 and [H? First, 
our systems are significantly larger: Refs, fl^ andfl^ used 
systems with up to 32^ and 128^ sites, respectively, while 
our systems have up to 1600^ sites. As the Ashkin-Teller 
model crosses over very slowly from the first-order tran¬ 
sition of the clean problem to the continuous transition 
of the disordered one, simulations of smaller systems are, 
perhaps, not sufficient to reach the asymptotic regime. 
Large systems are particularly important if the disorder 
strength is small. In fact, our own simulations for the 
weak dilution p = 0.1 show that, depending on e, the 
asymptotic regime may only be reached for L > 400. 
The random-bond system with Jh = 6/5 and J; = 4/5 
is especially weakly disordered; correspondingly, it does 
not reach the asymptotic regime even for L = 1600. 

This interplay between the disorder strength and the 
cross-over between first-order and continuous transitions 
is also borne out by the analysis of the correlation length 
exponent v. The asymptotic finite-size scaling (1141) of 
dimensionless quantities such as the Binder cumulant 
leads to an effective exponent r'eff > 1- This is what 
we have observed in all our systems except for the one 
with the weakest disorder, viz., the random-bond system 
with with Jh = 6/5 and Ji = 4/5. This supports the 
notion that the correlation length exponents reported in 
Refs. [H and [H may be effective exponents outside the 
asymptotic regime. 

We note, however, that a signihcant discrepancy be¬ 
tween our data and those reported in Ref. |T^ is manifest 
already in the clean phase diagram. Fig. [31 where system 
size effects should be less important. Our phase bound¬ 
ary agrees with the old data by Grest and Widon>i^ but 
disagrees with Ref. M 

As a byproduct, our simulations for e = 0 (where the 
Ashkin-Teller Hamiltonian is equivalent to three indepen¬ 
dent Ising models) also help to resolve the long-standing 
controversy about the critical behavior of the disor¬ 
dered two-dimensional Ising model. Our large-scale data 
for systems with up to 2240^ sites provide strong sup¬ 
port for the logarithmic-correction (strong-universality) 
scenario2^“— according to which the critical behavior is 
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characterized by the clean Ising exponents and universal 
logarithmic corrections. 

We now put our results in the general context of 
phase transitions of two-dimensional disordered sys¬ 
tems. Following the analytical results on the disor¬ 
dered two-dimensional Ising2^“— and Ashkin Teller- 
models, it was conjectured that all critical behavior in 
two-dimensional disordered systems belongs to the dis¬ 
ordered Ising universality class. This belief in super- 
universal critical behavior was further strengthened by 
early numerical results for disordered Ising Ashkin- 
Tellerf^^ and Pott o^^i^^ models as well as heuristic inter¬ 
face argumentsHowever, later simulations of the disor¬ 
dered g-state Potts model^i^ belied these expectations: 
They showed that the exponent jijv does depend on the 
value of q and generally differs from the Ising value of 
1/8. Recently, unexpectedly complex behavior was also 
found in the two-dimensional random-bond Blume-Capel 
model^^— an Ising-like spin-1 model with an additional 
single-ion anisotropy. 

Cardy’s renormalization group approacbii was gener¬ 
alized by Pujol^ from N coupled Ising models to N cou¬ 
pled g-state Potts models. For q = 2 (the Ising case), 
Pujol’s results agree with Cardy’s. For q > 2, however, 
he found the emerging critical behavior to be controlled 
by a nontrivial random fixed point. Testing these predic¬ 
tions numerically remains a task for the future. 

Finally, the quantum version of the Ashkin-Teller 
model has recently attracted considerable attention in 
connection with the question of how first-order quan¬ 
tum phase transitions react to disorder. Strong- 
disorder renormalization group calculations predict 
infinite-randomness critical points in different universal¬ 
ity classes, depending on the coupling strength 
Moreover, the two-color model is predicted to feature an 
unusual strong-disorder infinite-coupling phased. Our 
Monte Carlo method can be easily generalized from the 
two-dimensional classical case to the (1 -I- I)-dimensional 
quantum case. Some calculations along these lines are 
under way. 
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Appendix A: Short Monte Carlo runs and unbiased 
estimators 


Short Monte Carlo runs consisting of only a small 
number of measurements per sample introduce biases 
into some observables, at least if one employs the usual 
estimators. Consider, for example, the magnetic sus¬ 
ceptibility (of a single sample) which is related to the 
variance of the magnetization via x — with 

cr^ = (m^) — (m)^. In a Monte Carlo simulation, the 
variance cr^ is usually replaced by the estimator 
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where mi is the magnetization of an individual measure¬ 
ment and Um is their number. It is well known in statis¬ 
tics that underestimates the variance, even for uncor¬ 
related mi. This can be seen by evaluating the expecta¬ 
tion value of as 
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If the mi are correlated with a correlation time 
of T, the bias becomes even stronger: (sm) ~ 
(1 “ (1 + ^)l'nm) with A T. Other quantities de¬ 
fined as variances or covariances develop analogous bi¬ 
ases, including the specific heat C = (L^/r^)((e^) —(e)^). 

Note that these biases are not important in normal 
Monte Carlo simulations that consist of one (long) run of 
n„i measurements because the bias decays as while 
the statistical error decays as nm ■ The bias is thus 
much smaller than the statistical error and can be ne¬ 
glected. However, if the results of short runs are aver¬ 
aged over a large number of samples, this argument 
changes. The bias still decays as but the statisti¬ 
cal error and sample-to-sample fluctuations due to disor¬ 
der are suppressed by an additional factor Us It is 
thus clear that the bias cannot be neglected for a suffi¬ 
ciently large number of samples. (If the disorder-induced 
sample-to-sample fluctuations are weaker than the ther¬ 
modynamic fluctuations, this is expected when Us > rim- 
In the opposite case, for strong disorder fluctuations, the 
bias becomes important roughly when Us > n^.) 

How can one correct the bias due to short Monte Carlo 
runs? If the measurements were completely independent, 
one could simply multiply the usual estimator IMl) by 
n-m/inm ~ !)■ However, achieving full independence re¬ 
quires long time intervals between consecutive measure¬ 
ments which makes the simulations inefficient. We in¬ 
stead introduce modified, unbiased estimators. To this 
end, we split the Monte Carlo run of Um measurements 
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into two halls, each with nml‘2‘ measurements. We also 
perform a few extra Monte Carlo sweeps between the two 
halfs to ensure that they are independent of each other. 
The improved estimator of tr^ is then given by 
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(A3) 

Following the same steps as outlined in (IA2I) . it is 
straightforward to show that (s^) = This means 
is unbiased. An analogous unbiased estimator can be 
defined for the specific heat C. 


In Sec. IIII D[ we defined two magnetic Binder cumu- 
lants, ^av and ggi, as well as two correlation lengths, ^av 
and ^gi. The “average” versions ^av and ^av suffer from 
short-run biases similar to those discussed above while 
the “global” versions gg\ and ^gi are unbiased. In princi¬ 
ple, one could correct the biases in ^av and ^av by using 
improved estimators. However these would have a more 
complicated structure than (IA3I) to deal with the terms 
in the denominators of eqs. and ©• For simplicity, 
we have not done this. Instead we mostly rely on the 
unbiased observables (?gi and ^gi. (Interestingly, our nu¬ 
merical data suggest that g^v and ^av have significantly 
smaller biases than C and x-) 
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